#  Replication code for analysis for
#  "The Limited Value of a Second Opinions"
#  William Minozzi and Jonathan Woon
#  2019/06/06

# Setup environment ----
suppressPackageStartupMessages({
  library(data.table)
  library(cowplot)
  library(rstan)
  library(rstanarm)
  library(extrafont)
  library(Cairo)
  })
load("results/data.RData")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
loadfonts(quiet = TRUE)
set.seed(360222565)
theme_set(theme_minimal() + theme(
  strip.background = element_rect(fill = alpha("white", 0.0),
    color = alpha("white", 0.0)),
  axis.title = element_text(size = 10), axis.text = element_text(size = 8),
  axis.ticks.x = element_blank(), plot.title = element_text(hjust = 0.5),
  panel.border = element_rect(color = "black", size = .2, fill = NA),
  text = element_text(family = "Linux Libertine"),
  strip.text = element_text(size = 10), legend.position = "bottom",
  panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank()))

# Estimate equilibrium action models ----
ad <- copy(Action_Data)
ad[Single == 1, Low_Region := as.numeric(Target <= -80)]
ad[Aligned == 1, Low_Region := as.numeric(Target <= -80)]
ad[Opposed == 1, Low_Region := as.numeric(Target <= 20)]
equilibrium_action_models <- list(
  stan_lmer(Action ~ Target * Low_Region + (1 | Session) +
      (1 | ID),
    data = ad[Single == 1], seed = 487536605,
    prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
      5), prior_covariance = decov(), prior_PD = FALSE,
    algorithm = "sampling", adapt_delta = NULL, QR = FALSE),
  stan_lmer(Action ~ Target * Low_Region + (1 | Session) +
      (1 | ID),
    data = ad[Aligned == 1], seed = 1779394190,
    prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
      5), prior_covariance = decov(), prior_PD = FALSE,
    algorithm = "sampling", adapt_delta = NULL, QR = FALSE),
  stan_lmer(Action ~ Target * Low_Region + (1 | Session) +
      (1 | ID),
    data = ad[Opposed == 1], seed = 1120466369,
    prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
      5), prior_covariance = decov(), prior_PD = FALSE,
    algorithm = "sampling", adapt_delta = NULL, QR = FALSE))
save(equilibrium_action_models,
  file = "results/equilibrium_action_models.RData")

# Estimate payoff models ----
payoff_model <- stan_lmer( (320 - Target_Action_distance) / 20 ~ Treatment +
    (1 | Session) + (1 | ID), data = Action_Data, seed = 999378286,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
payoff_period_model <- stan_lmer((320 - Target_Action_distance) / 20 ~
    Treatment * I(Period > 12) +
    (1 | Session) + (1 | ID), data = Action_Data, seed = 1011586230,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
save(payoff_model, payoff_period_model, file = "results/payoff-models.RData")

# Estimate width model ----
md <- copy(Message_Data)
md[, Second_Half := as.numeric(Period > 12)]
width_model <- stan_lmer(Message_Range_Size ~
    Moderate_Aligned + Extremist_Aligned + Moderate_Opposed +
    Extremist_Opposed +
    (1 | Session) + (1 | ID), data = md, seed = 174319797,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
width_time_model <- stan_lmer(Message_Range_Size ~
    Moderate_Aligned * Second_Half +
    Extremist_Aligned * Second_Half +
    Moderate_Opposed * Second_Half +
    Extremist_Opposed * Second_Half +
    (1 | Session) + (1 | ID), data = md, seed = 1219229678,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
save(width_model, width_time_model, file = "results/width-model.RData")

# Estimate sequential message models ----
ad <- copy(Action_Data)
ad[Single == 1, Low_Region := as.numeric(Target <= -80)]
ad[Aligned == 1, Low_Region := as.numeric(Target <= -80)]
ad[Opposed == 1, Low_Region := as.numeric(Target <= 20)]
seq_message_aligned_model <- stan_lmer(Moderate_Message ~
    Extremist_Message + Target +
    (1 | Session) + (1 | ID), data = ad[Treatment == "Aligned"],
  seed = 1664380463,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
seq_message_opposed_model <- stan_lmer(Moderate_Message ~
    Extremist_Message + Target +
    (1 | Session) + (1 | ID), data = ad[Treatment == "Opposed"],
  seed = 1417379256,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
seq_message_joint_model <- stan_lmer(Moderate_Message ~
    Extremist_Message * Opposed  +
    Target * Opposed +
    (1 | Session) + (1 | ID), data = ad[Treatment != "Single"],
  seed = 263015834,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
seq_message_aligned_hi_model <- stan_lmer(Moderate_Message ~
    Extremist_Message + Target +
    (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned" & Low_Region == 0],
  seed = 1664380463,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
seq_message_opposed_lo_model <- stan_lmer(Moderate_Message ~
    Extremist_Message + Target +
    (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed" & Low_Region == 1],
  seed = 1417379256,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
seq_message_opposed_hi_model <- stan_lmer(Moderate_Message ~
    Extremist_Message + Target +
    (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed" & Low_Region == 0],
  seed = 1417379256,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
save(seq_message_aligned_model, seq_message_opposed_model,
  seq_message_joint_model,
  seq_message_aligned_hi_model,
  seq_message_opposed_lo_model,
  seq_message_opposed_hi_model,
  file = "results/seq-message-models.RData")

# Estimate message models ----
ad <- copy(Action_Data)
ad[, Second_Half := as.numeric(Period >= 13)]
ad[, Moderate_Message_Midpoint :=
    (Moderate_Message_Min + Moderate_Message_Max) / 2]
ad[, Extremist_Message_Midpoint :=
    (Extremist_Message_Min + Extremist_Message_Max) / 2]
message_midpoint_single_time_model <- stan_lmer(Moderate_Message_Midpoint ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Single"], seed = 1119756487,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_min_single_time_model <- stan_lmer(Moderate_Message_Min ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Single"], seed = 1997393497,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_max_single_time_model <- stan_lmer(Moderate_Message_Max ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Single"], seed = 18208788,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_midpoint_aligned_extremist_time_model <- stan_lmer(
  Extremist_Message_Midpoint ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned"], seed = 1876123558,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_min_aligned_extremist_time_model <- stan_lmer(Extremist_Message_Min ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned"], seed = 584754035,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_max_aligned_extremist_time_model <- stan_lmer(Extremist_Message_Max ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned"], seed = 744110765,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_midpoint_opposed_extremist_time_model <- stan_lmer(
  Extremist_Message_Midpoint ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed"], seed = 501131098,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_min_opposed_extremist_time_model <- stan_lmer(Extremist_Message_Min ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed"], seed = 124616192,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_max_opposed_extremist_time_model <- stan_lmer(Extremist_Message_Max ~
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed"], seed = 609445167,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_midpoint_aligned_moderate_time_model <- stan_lmer(
  Moderate_Message_Midpoint ~
    Extremist_Message * Opposed * Second_Half  +
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned"], seed = 629295434,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_min_aligned_moderate_time_model <- stan_lmer(Moderate_Message_Min ~
    Extremist_Message * Opposed * Second_Half  +
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned"], seed = 629301878,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_max_aligned_moderate_time_model <- stan_lmer(Moderate_Message_Max ~
    Extremist_Message * Opposed * Second_Half  +
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Aligned"], seed = 264158998,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_midpoint_opposed_moderate_time_model <- stan_lmer(
  Moderate_Message_Midpoint ~
    Extremist_Message * Opposed * Second_Half  +
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed"], seed = 1108614329,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_min_opposed_moderate_time_model <- stan_lmer(Moderate_Message_Min ~
    Extremist_Message * Opposed * Second_Half  +
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed"], seed = 1928230970,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
message_max_opposed_moderate_time_model <- stan_lmer(Moderate_Message_Max ~
    Extremist_Message * Opposed * Second_Half  +
    Target * Second_Half + (1 | Session) + (1 | ID),
  data = ad[Treatment == "Opposed"], seed = 1710912810,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
save(
  message_midpoint_single_time_model,
  message_min_single_time_model,
  message_max_single_time_model,
  message_midpoint_aligned_extremist_time_model,
  message_min_aligned_extremist_time_model,
  message_max_aligned_extremist_time_model,
  message_midpoint_aligned_moderate_time_model,
  message_min_aligned_moderate_time_model,
  message_max_aligned_moderate_time_model,
  message_midpoint_opposed_extremist_time_model,
  message_min_opposed_extremist_time_model,
  message_max_opposed_extremist_time_model,
  message_midpoint_opposed_moderate_time_model,
  message_min_opposed_moderate_time_model,
  message_max_opposed_moderate_time_model,
  file = "results/message_models.RData")

# Esimate Action model ----
action_model <- stan_lmer(Action ~ Moderate_Message + Extremist_Message +
    Moderate_Message : Aligned + Moderate_Message : Opposed +
    Extremist_Message : Opposed + Aligned + Opposed +
    (1 | Session) + (1 | ID), data = Action_Data, seed = 440033593,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
action_extremist_model <- stan_lmer(Action ~ Extremist_Message +
    Extremist_Message : Opposed + Opposed + (1 | Session) + (1 | ID),
  data = Action_Data[Treatment != "Single"], seed = 174319797,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
action_moderate_model <- stan_lmer(Action ~ Moderate_Message +
    Moderate_Message : Aligned + Moderate_Message : Opposed +
    Aligned + Opposed + (1 | Session) + (1 | ID),
  data = Action_Data, seed = 1664380463,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
action_first_half_model <- stan_lmer(Action ~
    Moderate_Message + Extremist_Message + Moderate_Message : Aligned +
    Moderate_Message : Opposed + Extremist_Message : Opposed +
    Aligned + Opposed + (1 | Session) + (1 | ID),
  data = Action_Data[Period <= 12],
  seed = 1417379256,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
action_second_half_model <- stan_lmer(Action ~
    Moderate_Message + Extremist_Message + Moderate_Message : Aligned +
    Moderate_Message : Opposed + Extremist_Message : Opposed +
    Aligned + Opposed + (1 | Session) + (1 | ID),
  data = Action_Data[Period >= 13],
  seed = 1149454341,
  prior = normal(), prior_intercept = normal(), prior_aux = cauchy(0,
    5), prior_covariance = decov(), prior_PD = FALSE,
  algorithm = "sampling", adapt_delta = NULL, QR = FALSE)
save(action_model, action_extremist_model, action_moderate_model,
  action_first_half_model, action_second_half_model,
  file = "results/action-model.RData")

# Estimate SF messages model ----
md <- copy(Message_Data)
setorder(md, ID)
md[, Message_Midpoint := (Message_Min + Message_Max) / 2]
md[Message_Midpoint == 0, Max_N_Messages := 201]
md[Message_Midpoint > 0, Max_N_Messages := 1 + 2 * (100 - Message_Midpoint)]
md[Message_Midpoint < 0, Max_N_Messages := 1 + 2 * (100 + Message_Midpoint)]
md[, Message_Range_Size := Message_Max - Message_Min]
md[Message_Midpoint == 0, Max_Message_Range_Size := 200]
md[Message_Midpoint > 0,
  Max_Message_Range_Size := 2 * (100 - Message_Midpoint)]
md[Message_Midpoint < 0,
  Max_Message_Range_Size := 2 * (100 + Message_Midpoint)]
md[, id := as.numeric(factor(ID))]
md_1 <- md[!Message_Midpoint %in% c(100, -100) & N_Messages != Max_N_Messages]
md_2 <- md[!Message_Midpoint %in% c(100, -100) & N_Messages == Max_N_Messages]
md_3 <- md[Message_Midpoint == 100]
md_4 <- md[Message_Midpoint == -100]
dl_messages <- list(
  N_id = max(md$id),
  N_period  = max(md$Period),
  N_1 = nrow(md_1),
  id_1 = md_1$id,
  period_1 = md_1$Period,
  Target_1 = md_1$Target / 100,
  Message_Midpoint_1 = md_1$Message_Midpoint / 100,
  log_Message_Range_Size_1 = log1p(md_1$Message_Range_Size / 100),
  N_2 = nrow(md_2),
  id_2 = md_2$id,
  period_2 = md_2$Period,
  Target_2 = md_2$Target / 100,
  Message_Midpoint_2 = md_2$Message_Midpoint / 100,
  log_Max_Message_Range_Size_2 = log1p(md_2$Max_Message_Range_Size / 100),
  N_3 = nrow(md_3),
  id_3 = md_3$id,
  period_3 = md_3$Period,
  Target_3 = md_3$Target / 100,
  U = 1,
  N_4 = nrow(md_4),
  id_4 = md_4$id,
  period_4 = md_4$Period,
  Target_4 = md_4$Target / 100,
  L = -1)
sm_messages <- stan_model(file = "stan/messages.stan")
sf_messages <- sampling(sm_messages,
  chains = 4, iter = 2000, warmup = 1000,
  data = dl_messages, seed = 1619334482)
save(sf_messages, file = "results/sf_messages.RData")
